Temperature Dependence of Facet Ridges in Crystal Surfaces. 
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The equilibrium crystal shape of a body-centered solid-on-solid (BCSOS) model on a honeycomb 
lattice is studied numerically. We focus on the facet ridge endpoints (FRE). These points are 
equivalent to one dimensional KPZ-type growth in the exactly soluble square lattice BCSOS model. 
In our more general context the transfer matrix is not stochastic at the FRE points, and a more 
complex structure develops. We observe ridge lines sticking into the rough phase where the surface 
orientation jumps inside the rounded part of the crystal. Moreover, the rough-to-faceted edges 
become first-order with a jump in surface orientation, between the FRE point and Pokrovsky- 
Talapov (PT) type critical endpoints. The latter display anisotropic scaling with exponent 2 = 3 
instead of familiar PT value z = 2. 

PACS number(s): 68.35 Bs, 64.60 Fr, 64.60 Ht, 68.35 Rh 



I. INTRODUCTION 

Equilibrium crystal shapes typically consist of various 
flat facets connected by rounded surfaces. The way this 
structure changes with temperature has been studied for 
many years |l]-|| . Most aspects are well established. For 
example, rounded parts are associated with crystal ori- 
entations where the surface is rough, and flat facets can 
disappear from the equilibrium shape at the roughening 
temperature of those specific facets. Roughening transi- 
tions belong to the Kosterlitz-Thouless (KT) universality 
class |lC]] . Many details follow from the scaling properties 
of the KT transition, e.g., that the facet diameter van- 
ishes exponentially with reduced temperature, and that 
the curvature of the surface is universal just above the 
KT transition at that facet orientation. 

Flat facets typically connect smoothly to rounded, 
rough surface areas. The surface orientation in the 
rounded part varies continuously and connects without 
a jump in angle onto the flat facet. Such phase bound- 
aries are described by Pokrovsky-Talapov (PT) transi- 
tions jll] The angle difference vanishes with a | 
power. The most salient feature of PT transitions is their 
anisotropic scaling. The rounded surface on the rough 
side of the PT transition can be interpreted as a stepped 
surface from the perspective of the flat facet. The char- 
acteristic lengths along, £m, and perpendicular to these 
steps, scale as £m ~ with z = 2. 

Two dimensional equilibrium statistical mechanical 
systems are mathematically equivalent to one dimen- 
sional quantum systems at zero temperature. The trans- 
fer matrix formalism of the first is the path integral 
representation of the second. PT transitions are equiv- 
alent to metal-insulator transitions in one dimensional 
fermion systems. The steps represent the world lines of 
fermions, and the lattice direction parallel to the steps 
plays the role of time. In the rough phase (the metal 
phase) the fermions behave relativistically, i.e., time, 
scales compared to space, with dynamic exponent 
z = 1, £ii ~ At the PT metal-insulator transition the 



fermions behave non-relativistically, with dynamic expo- 
nent z = 2 instead of z = 1 HQ. 

The rough-to-facet interface can also be first-order, 
with a discontinuity in the surface orientation angle in- 
stead of the PT-type smooth connection. The critical 
endpoint between such first-order and PT-type segments 
provides an example of a quantum field theory with dy- 
namic exponent z = 3. To our knowledge an explicit 
example of such a PT critical endpoint (PTE) has not 
been reported before. We will do so in this paper. 

Our current interest in these issues stems from the 
discovery that the scaling properties at facet ridge end- 
points (FRE) are related to KPZ-type growth [[§. FRE 
points are the endpoints of ridges where two facets con- 
nect directly. The ridge splits into two PT or two first- 
order transitions with a rounded, rough surface in be- 
tween. The exactly soluble body-centered solid-on-solid 
(BCSOS) model on a square lattice (with only next near- 
est neighbour interactions) contains such a FRE point 
||. Moreover, the transfer matrix of this model becomes 
stochastic at the FRE point and maps exactly onto the 
master equation time evolution operator for the so-called 
brick model in one dimension lower |^, p^ , |l5| ] . The latter 
describes a growing ID crystal. The ID growing inter- 
face is rough, and its scaling properties belong to the 
KPZ universality class. The _L-direction represents the 
spatial direction and the ||-direction represents time, just 
as in mapping of the PT transtion to ID fermions men- 
tioned above. The dynamic exponent is equal to z = |. 
The characteristic time needed to reach the stationary 
growing state diverges as fn ~ N z , with N the size of 
the ID lattice. 

This translates into non-trivial scaling properties of 
KPZ-type FRE points When approached from the 

rounded phase in the direction parallel to the facet ridge, 
see Fig. [j], the surface orientation jumps. Moreover, at 
the rough side the characteristic length scale in the di- 
rection along the facet ridge, scales as £|| ~ with 

z = ^ 
a 2 . 

This exact mathematical mapping between the scaling 
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properties of equilibrium phase transitions to those of dy- 
namic non-equilibrium transitions in one lower dimension 
is intriguing. But it leaves us with an important ques- 
tion. The above BCSOS model satisfies a special sym- 
metry at its FRE point. The transfer matrix is "acciden- 
tally" stochastic. It has full rotational symmetry in the 
spin-i operators (the step variables). This implies, e.g., 
a complete degeneracy of the free energy with respect to 
all surface orientations. Such a symmetry is too special. 
FRE points arc a more general phenomenon, not limited 
to KPZ-type stochastic transfer matrices. We need to in- 
vestigate whether the KPZ-type scaling exponent z = | 
is generic. Maybe stochasticity is restored at large length 
scales, asymptotically close to the FRE point, with non- 
stochasticity an irrelevant operator. More likely, the FRE 
point changes its character. In this paper we present such 
an investigation, a numerical study of a model with a 
complex crystal shape. 

We define this model in section 2. It originates 
from our earlier study of possible spontaneous low-angle 
faceting in FCC(lll)-type facets This model has 

the required property that its crystal shape is complex 
and that accidental symmetries are unlikely. Section 3 
contains an overview of the phase diagram and crystal 
shapes. In section 4 we describe our numerical method. 

The crystal structure is actually much richer than we 
anticipated. Entire facets vanish as functions of the cou- 
pling constants. FRE points become triple points. We 
also find a more exotic feature, like a first-order (FOR) 
line inside the rough phase emerging from the FRE point 
(section 5). At this line the surface orientation angle in- 
side the rough phase jumps. In their paper about the 
exact solution of the FRE point in the square lattice BC- 
SOS model, Bukman and Shore || speculate about the 
possible generalizations of Fig. |l|. The FOR line structure 
corresponds to their so-called ridge scenario. However, 
they doubted that it would be realized. In our model, 
FOR lines are clearly present although they remain very 
short. 

In a nearby section of the phase diagram, the FOR lines 
disappear. The rough-to-faceted ridge becomes sharp 
with a jump in the surface orientation angle. These first- 
order lines terminate further away from the FRE point 
in PTE points and the edge continues as a second order 
PT line. In section 6 we study the scaling properties at 
these PTE points. They display anisotropic scaling with 
exponent 2 = 3. The conclusions are presented in section 
7. 



II. THE MODEL 

Our model emerges naturally at a coarse grained level 
in a study of spontaneous low angle faceting in FCC(lll) 
facets. This will be described elsewhere [fig . For clar- 
ity it is more appropriate to introduce the model from 
a different perspective. Consider the (001) facet of an 



HCP-type crystal, with ABABA-type stacking of trian- 
gular slices. The solid-on-solid (SOS) model description 
of such surfaces leads to a BCSOS model on a honey- 
comb lattice. The surface heights on the A-type triangu- 
lar sublattice are even integers, Ha — 0, ±2, ±4, . . ., and 
those on the f?-type triangular sublattice are odd inte- 
gers, Kb = ±1,±3,±5, . . . The simplest Hamiltonian for 
such a surface is 

(id) 

with the summation over next nearest neighbour 

columns. 

We break the symmetry between the A- and -B-type 
columns, such that two distinct (001) facets can be cut 
from this crystal, those with A's or those with B's on top. 
A and B could be different type of molecules, but that is 
not what we have in mind. The interactions between the 
next nearest neighbour columns, La and Lb, would then 
be different. Instead, assume that A and B are identical 
types of atoms, but that for some reason the bulk struc- 
ture is reconstructed such that the A's move downward 
from the positions exactly in-between the B layers. In 
that case the bonds between the A atoms with the B's 
in the layer beneath them are stronger than with the B's 
in the layer above them. We can model this by adding a 
nearest neighbour coupling in the BCSOS model 

H = K ^(hi-h^ + ^L^ihi-hj) 2 (2) 

<id> (id) 

where, in the sum over nearest neighbors, i (j) always 
refers to A (B) sublattice. 

An alternative more generic scenario is the one real- 
ized in Naphthalene j|] (for a different crystal symmetry 
however). It involves only a surface reconstruction, no 
bulk reconstruction. Consider a crystal where the A's 
and B's are identical molecules, but with an orientation 
degree of freedom pointing in two different directions. 
Inside the bulk these two orientations can be equiva- 
lent by (a slide-mirror-plane type) symmetry, but the 
presence of the surface typically breaks that equivalence. 
The surface molecules rotate due to surface relaxation. 
Typically they do so differently with the A's or the B's 
on top. This breaks the A B symmetry locally close 
to the surface and creates an energy difference, K ^ 0. 
The next nearest neighbour interactions remain virtually 
equal, La — Lb = L, if the orientational aspects of the 
interactions decay rapidly with distance. 

We introduce tilt energies Ei, i — 1,2, 3, to study the 
crystal structure at and nearby the central (001) facet 
orientation. The functional dependence of the free en- 
ergy with respect to these Ei yields the crystal shape by 
means of the Wulff construction || . 

n= \LY, {id) {hi-h j ) 2 + 

(K + E 1 )J2 <ij>i (h i -h j ) + 
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(K+E 2 )j: <i>j>2 (h i -h j ) + 



(3) 



The summations run over nearest neighbour columns in 
the three different directions, see Fig. [| The three tilt 
fields Ei are related as 



E± — E x 

~2 y 
2 



E2 — -^E qj 
E 3 = — — E,. 



— Erf 
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(4) 



E x and E y are conjugate to the tilt angles in the x- and 
y-directions, i.e., parallel and perpendicular to one of the 
three main axes of the honeycomb lattice. 

We simplify the model further by not allowing step ex- 
citations to touch each other. Fig. ^ shows all the remain- 
ing local configurations in this 10-vertex model. The ar- 
rows denote whether nearest neighbour AB columns dif- 
fer in height by +1 or — 1. The thin (fat) lines denote that 
the K bond is in the low (high) energy state. The Boltz- 
mann weight contributions are shown in Fig. [| below each 
vertex configuration, with z c = e~ 2K ~ 2L , z s — e~ 



-2K-4L 



= e- B */ 4 ,and. 



-V3E y /A^ rpk e C0U pij n g cons tants 



are made dimensionless by absorbing the factor ksT. 



III. THE PHASE DIAGRAM 

The phase diagram of our model at zero tilt fields 
E x = E y = 0, Eq.(|), is shown in Fig. |. The (001) facet 
is stable at all temperatures inside the quadrant where 
K and L are both positive. The step excitations do not 
create sufficient entropy to roughen this surface even at 
T = 00 (at point K = L = 0). The roughening line lies 
inside the other three quadrants. We determine its loca- 
tion in the conventional manner for transfer matrix finite 
size scaling calculations [jl7|. The step free energy 77 in 
semi-infinite strips of width TV is equal to the difference 
in free energy for periodic and stepped boundary condi- 
tions in the finite direction. The KT roughening line is 
obtained by extrapolating the lines where Nr) = jir for 
increasing strip widths. 

At low temperatures for L < the surface sponta- 
neously facets into three coexisting orientations. We refer 
to these as straight-step (SS) facets, shown in Fig. [s](a)to 
distinguish them from the three zig-zag (ZZ) structures, 
shown in Fig. ||(b), that appear elsewhere in the phase 
diagram at non-zero tilt fields Ei. The transition from 
the rough phase into this phase with three coexisting SS 
facets is first-order. The surface orientation angle jumps. 
The location of this transition line is determined by the 
methods described in section 4. Its location is not triv- 
ial, unlike the square lattice BCSOS model || where it is 
simply the meeting point of several Ei ^ PT transition 
lines at Ei = 0. 



An additional ground state appears in the quadrant 
with K < and L > 0, the triangles-state shown in 
Fig. ||(c). This state represents a highly degenerate 
surface reconstruction. Unfortunately this ground state 
is unstable against thermal fluctuations immediately at 
T > in this particular model. The surface reconstruc- 
tion and associated preroughening phenomena associated 
with this type of ground state could have been quite in- 
teresting. 

We determine the crystal shape at each value of K and 
L by introducing the tilt fields Ei, see Eq. (||). The tech- 
nical details of these calculations are described in section 
4. In this section we summarize the crystal shape de- 
velopment with temperature. We identified at least 5 
distinct crystal structures, labeled I through V in Fig. ||, 
which shows at which values of K and L they are realized, 
as an overlay of the zero-tilt phase diagram Fig. [|. 

In region I, located at L > 0, the crystal structure 
consists of a non-tilted flat facet surrounded by a rough, 
rounded central region, and fully tilted facets of both 
types, SS and ZZ, at the edges (each repeated three 
times). This structure is shown in Fig. All facets are 
separated from the rough region by PT lines. Their exact 
location can be calculated easily for the SS facets. Com- 
pare the free energy density of the faceted phase with 
that of a state with one line defect, e.g., as shown in 
Fig. The PT transition occurs when the free energy of 
the faceted phase becomes equal to the free energy with 
the defect. The calculation is straightforward, but leads 
to a somewhat complex formula. 

The PT lines in region I never touch or join. There 
are no FRE points between the SS and ZZ facets. The 
following analysis explains why. Suppose facet ridges be- 
tween the ZZ and SS facets exist. Both faceted structures 
are frozen, without any fluctuations (in our model) and 
therefore the free energy in both is equal to the energy, 
E zz = 2K + 2L- \E X and E ss = 2K + 4L - \E X - 
\\fiEy. At a facet ridge these must be equal. This 
happens at E y = -^L. However, the ZZ and SS facets 
roughen before this line is reached. Consider the zero 
temperature limit of the Ezz = Ess line. Line defects 
such as the one shown in Fig. || have the same energy 
as the two faceted configurations. The facets are un- 
stable with respect to these defects, because at non-zero 
temperatures they gain entropy while the facets remain 
frozen. Two PT lines emerge from zero temperature and 
preempt the facet ridge. 

As L decreases, the ZZ facets shrink and disappear at 
L = 0. At L < 0, the only remaining fully tilted facets 
are the three SS ones. Along L = 0, they are still com- 
pletely separated from each other, with rough, rounded 
orientations in between them, and PT transitions as bor- 
ders. For L < 0, facet ridges appear between the three SS 
facets. These facet ridges are capped by FRE points. The 
local structure at these FRE points changes with temper- 
ature in an exotic manner as shown in Fig. ^. The global 
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crystal shape structure is less interesting. Initially it con- 
tains the central flat facet but this shrinks and vanishes 
across the roughening line. 

Fig. ^|(a) shows the local structure around the FRE 
point in region II. Two PT lines emerge from the FRE 
point and form the borders between the SS facets and 
the rough phase. Moreover, a FOR line emerges from the 
FRE point inside the rough phase. At this line the crys- 
tal orientation angle jumps. The skipped unstable sur- 
face orientations are associated with local ZZ facet type 
zig-zag step configurations. Those are unfavourable for 
L < and large K. The FOR line is definitely resolved 
in our numerical data, but remains extremely short (sec- 
tion 5). At its longest extension, it reaches only about 
two percent of the distance from the FRE point to the 
center. In some parts of region II it is so short that it 
can only barely be distinguished numerically, if at all. 

The FOR lines disappear at the boundary between re- 
gions II and III. This boundary looks sharp in Fig. 
but in reality it must be a smooth type of crossover. We 
discuss this in section 6. 

In region III, see Fig. ^, the FOR line vanished. 
The border between the SS facets and the rough phase 
changed as well. Close to the FRE point the SS-rough 
boundaries are now first-order. They continue as PT lines 
beyond PTE- type critical endpoints. Fig. ||(b) shows this 
local structure. 

The first-order segments become longer towards region 
IV. At the III- IV phase boundary, two PTE points from 
opposite FRE sides meet, and in region IV, the entire SS- 
rough phase boundary is first-order. In the mean time, 
the central rough region has been shrinking, and it dis- 
appears completely at the IV-V phase boundary. 



IV. NUMERICAL METHODS 

We determine the crystal shape by means of the trans- 
fer matrix method. This means that we calculate the 
exact free energies for semi-infinite lattices in terms of 
the largest and nearby eigenvalues of the transfer ma- 
trix. This is a standard method. Therefore we only need 
to comment here on the details concerning the surface 
tilt fields Ei, the amount of surface tilts Qi they cre- 
ate, and the Wulff construction (Legendre transforma- 
tion) between these two types of quantities. 

There are two natural choices for the directions in 
which the lattice is infinite. T can be aligned with the 
y-axis or the x-axis, see Fig. ||. We refer to these as T x 
and T y type transfer matrices. In T x , the infinite lattice 
direction coincides with the y-axis and the finite direc- 
tion with the x-axis. In T y it is the other way around. 
Notice that for each choice there are three rotationally 
equivalent ones. For T x , we are able to find the largest 
eigenvalues for 4 < N < 18, where N is the number of 
vertices in the time slice. For T y , we can calculate the 
largest eigenvalues for 2 < N < 10. 



Consider the free energy as function of tilt angles 
f{Qxi Qy)- This is the free energy of a facet at a given ori- 
entation (Q x ,Q y ). Some orientation ranges are thermo- 
dynamically unstable. These angles are not represented 
in the crystal shape. They are skipped, and result in 
sharp edges in the surface. To obtain the crystal shape 
one needs to minimize f{Q x ,Q y ) under the constraint 
that the volume of the crystal (the amount of matter) 
is conserved. This leads to the famous Wulff construc- 
tion [Q which is in essence a geometric construction for 
the Legendre transformation of f(Q x , Qy)- The tilt fields 
Ei are the conjugate variables to the tilt angles Qi. The 
crystal shape is a direct representation of the free energy 
function f(E x ,E y ) at a certain temperature (at specific 
values of K and L) |J. The shape functions Q X (E X , E y ) 
and Q y (E x , E y ) tell us which surface orientations are rep- 
resented in the crystal. 

The transfer matrix provides us with a mixed version of 
this. It leads to f(Q±, En) and Qu (Qx, En). We control 
the tilt angle in the finite lattice direction, Q±, and the 
tilt field in the infinite direction, Eu . Qy sets itself. Q± is 
set by the boundary condition in the finite lattice direc- 
tion, h{x\\ , x± + N) — h(xu,x±) + a. The total tilt of the 
surface can take the values a = 0, ±2, ... up to TV when 
the transfer matrix is aligned with the y-axis. This means 
that an average slope Q x = ^/3a/2N is frozen into the 
surface in the x direction. In the opposite set-up, with 
T y , where the lattice is infinite in the x direction, a tilt 
Q y = a/2N is frozen into the surface in the y-direction. 

Consider the transfer matrix at zero tilt field E± = 
and a specific value of Eu . We determine the largest 
eigenvalue, Ao of T for each value of Q±. The free en- 
ergy density is equal to f(E\\,0,Qx) = - log(A )/iV. 
Fig. [l0|(a) shows a typical example, at K = —0.24, 
L = 0.29, and Eu = 0. Only the Q± > side is shown, 
because the curve is mirror symmetric. The discrete set 
of possible tilt angles Q± increases with the strip width 
N, and we need to perform a finite size scaling (FSS) 
analysis to obtain the infinite- by-infinite lattice result. 
Fortunately these FSS corrections are reasonably small, 
and typically converge smoothly. 

The tilt angle Q± that is realized is the one that min- 
imizes the free energy. This corresponds to the low- 
est point in Fig. |lC](a). There is no need to repeat 
the calculation for other values E± ^ 0. They are 
related by f(E {] ,E ± ,Q ± ) = f(E {] , 0, Q ± ) + E ± Q ± be- 
cause Q± is conserved by the transfer matrix. In other 
words, the crystal shape f(E±,Eu) and the shape func- 
tion Q±(E±, E\\) follow from the Legendre transform of 
/(i?H , with respect to Q± 

It is easy to see from graphs like Fig. |l^(a) how Q± 
varies with Ej_ . Simply rotate the page slowly and watch 
the minimum shift, and sometimes skip certain orienta- 
tions. In Fig. |l^(a) such skips do not take place. We 
follow the iS^-axis at E\\ = E y = in Fig. [?] along which 
there are no orientational discontinuities. Notice that 
f{Eu = 0, Q±) develops a cusp at Q x = with system 
size. This means that the central (001) facet is stable at 
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small E±. The (OOl)facet boundary is a PT transition 
because the cusp is convex at Q± = and no orientations 
are being skipped for increasing E± (rotate the page). 
The KT roughening transition takes place at values K 
and L where the cusp disappears. In the rough phase, 
/(-Ell, Q±) has a quadratic minimum at Q± = 0. Its cur- 
vature represents the roundness of the crystal shape at 
that orientation. 

Next, consider the stability of the SS fully tilted facet. 
At maximum tilt angle the free energy f(E»,Qx) has 
a definite slope. This means that the SS facet is sta- 
ble for a specific range of E± = E x . The absence of 
such a slope would imply that the SS facet is thermody- 
namically unstable, and would be absent in the crystal 
shape. The SS-rough phase boundary is a PT transition 
because f(E\\,Q±) is convex at Q± near its maximum 
and no values of Qj_ are being skipped by the Legendre 
transformation. 

The above analysis is still incomplete because we did 
not check how Q» is shifting simultaneously. The func- 
tion Q II (.En, Q±) is contained in the right and left eigen- 
vector of the transfer matrix assocated with the largest 
eigenvalue. In the interpretation of T as a quantum me- 
chanical time evolution operator, Qn is a velocity type 
off-diagonal expectation value. 

Fig. |l(](b) shows Q|| {E\\ , Q±) = Q y {E y , Q x ) at our ex- 
ample point. It shows that the orientation of the step 
excitations in the surface rotates as expected, from verti- 
cal just outside the Q± — Q\\ = flat facet to the correct 
slope just outside the SS fully tilted facet. 

In the above example, we used the T x transfer matrix 
set-up to investigate the crystal shape in the horizontal 
direction. In that set-up Q x takes only discrete values, 
Qx = V3a/2N. The opposite set-up with T y would have 
been better, since then Q x would have been continuous 
and Q y would have been rational. Fig. |ll| illustrates this. 
We still use T x but now determine the crystal structure 
along the E y axis. The free energy curve lacks mirror 
symmetry. It shows the SS facet at the Em < side and 
the ZZ facet at the Ey > side. It also shows the central 
facet around E^ = 0. The facet edges are smooth since 
they are PT-like. 

From the above example it must be clear how we deter- 
mined the development of the crystal shape with K and 
L, described in section 3. It required numerous graphs 
of the type shown in Fig. [II] and Fig. |ll|. The results 
are quantitively represented in Fig. ||, Fig. (?], and Fig. |j. 
The latter two show the crystal shapes at actual values 
of K and L as specified in the figure captions. In the fol- 
lowing two sections we limit the discussion to the most 
interesting features of the crystal shapes: the FOR lines, 
and the PTE points (see FigM. 



V. REGION II: THE FOR LINES 



In region II the crystal shape contains only SS-type, 
fully-tilted facets, with facet ridges in between them, 
capped by FRE points. The only structural difference 
with the KPZ-type FRE points in the square lattice BC- 
SOS model is the presence of the FOR lines sticking into 
the rough phase, compare Fig. ^(a) with Fig. [|. 

Such first-order ridges inside the rough phase have 
never been realized theoretically or experimentally be- 
fore. Unfortunately, they remain very short in our model 
and are difficult to resolve numerically. Fig. ^(a) is drawn 
on scale for point K = 1.43 and L — —0.63 where the 
FOR line is at its longest, but remains very short com- 
pared to the distance of the FRE point to the center, 
AE/E ~ 0.02. 

In the following we provide numerical evidence that 
FOR lines in our model are for real and not a fancy, 
transient, finite size scaling phenomenon. This is an im- 
portant issue. Suppose the FOR lines vanish in the ther- 
modynamic limit. Region II then represent a realization 
of the same type of FRE points as the ones in the square 
lattice BCSOS model, except that the transfer matrix is 
not stochastic anymore. We determined the anistropic 
scaling exponent z from the leading, Ao, and next lead- 
ing, Ai, eigenvalues of the transfer matrix in the sector 
with zero tilt Q± = 0. The massgap m = log(Ao/Ai) 
represents the inverse of the time like correlation length, 
£|| ~ m _1 ~ N z . We do this at our best estimates for the 
endpoint of the FOR line at several values of K and L, 
and find z ~ 1.6 ± 0.1 (see Fig. 12). This is close to the 
KPZ exponent z = 3/2. The assumption that the FOR 
lines vanish in the thermodynamic limit would confirm 
that KPZ scaling is valid beyond the stochastic subspace. 
However, in that interpretation, it remains unclear how 
far we are from the stochastic subspace. Crossover cor- 
rections to scaling might mask the "true" non-stochastic 
exponent z. 

The existence of the FOR lines is much more likely, 
however. Fig. 13 shows the variation of the free energy 
f(Q x ) at K 



1.43 and L = —0.63 for various E y close 
to the FRE point in Fig. ||(a). We need to use the T x 
set-up, the one with the infinite lattice direction aligned 
with the FOR line. The rational values of Q y in the T y 
set-up skip across the FOR line and make it invisible al- 
together (since it is so short). In the T x set-up the Q x are 
rational. In Fig. |l3| we show only the data for N = 18, 
otherwise the graph becomes too crowded. 

Fig.[l4| shows the corresponding Q y (Q x ) behaviour at 
one value of E y (the curves for the other values of E y 
along the FOR line lie almost on top of this curve) . Well 
inside the rough phase (below the FOR line in Fig. [|(a)), 
the f(Q x ) curves are convex. The minimum is at Q x = 0, 
and it shifts smoothly with E x (rotate the page). More- 
over, the curve remains convex all the way into the wings, 
implying PT-type boundaries with the SS facets. 

In the opposite limit, well above the FRE point, the 
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curves are concave. The two SS facets are coexisting at 
E x = 0, the entire f(Q x ) curve except for its boundary 
values in the wings represent thermodynamically unsta- 
ble orientations. 

The intermediate behaviour of the f(Q x ) curve deter- 
mines the crystal shape close the FRE point. In the 
square lattice BCSOS model it flips at once from com- 
pletely convex to completely concave (like changing the 
sign of a in a polynomial like f(Q x ) — a Q x )- At that 
FRE point the transfer matrix is stochastic, implying 
that all surface orientations have the same free energy. In 
general, f(Q x ) has more structure. There are many pos- 
sible structures, but we observe in our model the simplest 
ones. In Fig. |l3|, the curvature of the central part of the 
f(Qx) curve changes before its SS facet wings come down. 
This creates the W-type shapes, and therefore gives rise 
to the FOR line. All of this happens in a very narrow E y 
interval. In Fig. [l5]we show the location of E y where the 
curvature at the central part becomes zero, as function 
of N for K = 1.43 and L — —0.63. To be more precise, 
we show the location of Xq(Q x ) = Ao(Q x + V3/2N) for 
Q x = 0, Q x = VS/2N, and Q x = 2^/S/2N. The curves 
pass the FRE point at small N, but then bend backward. 
Still, they do not converge to the FRE point (its location 
is known exactly and shown in the figure for reference). 
The curve for Q x = is straightening out for system sizes 
12 < N < 18 and converges to a point below the FRE 
point at E y — 0.910 ±0.005. We conclude from this that 
the FOR lines must be for real. 



VI. REGION III: SCALING AT PTE POINTS 

In region III the f(Q x ) curves near the FRE point are 
M shaped (inverted W's). This is shown in Fig.|l6| for 
K = 0.37 and L = —0.26. Fig.[l7] shows the correspond- 
ing Q y (Q x ) plots. The SS wings of the f(Q x ) curve bend 
down and cross the central minimum before the curva- 
ture at Q x = has a change to change sign. This be- 
haviour implies that the f{Q x ) curves become concave 
in the wings, and therefore that the PT transitions be- 
come first-order, see Fig. 0(b). Contrast this with the W 
shapes in region II, Fig. O , where the curvature of the 
central part of f(Q x ) changes sign before its wings come 
down, and thus creates the FOR line. 

The phase boundary between regions II and III, can 
not be sharp. It looks abrupt in Fig. ^, but that is an ar- 
tifact of the shortness of the FOR lines. This changeover 
could take many forms. The FOR lines require the devel- 
opment of a concave segment in the center of f(Q x , E y ) 
while the PTE points require the development of concave 
segments in the wings. These two phenomena could be 
independent, and that would result in both a FOR line 
and first-order SS-rough segments, as shown in Fig.fi8|(a). 
The FRE points become first-order as well. This means 
that subtle developments in the central part will be 
skipped and not be expressed in the crystal shape if they 



take place after the abrupt crossover from the central 
minimum to the wing minima. Our FOR lines are very 
short and move rapidly into the thermodynamically un- 
stable region to vanish from the crystal shape. But this 
is only one of the possible scenarios. 

We can not resolve accurately how the structure 
evolves from II to III. It happens too fast and the FOR 
lines are too short. We find some evidence that two con- 
cave bubbles move rapidly along the f(Q x ,E y ) curves 
from the center to the wings in opposite directions. In 
that case the FOR line splits into two first-order lines, 
as shown in Fig. [l8| Each line is still inside the rough 
phase, and each merges with the PT lines to form the 
first-order SS-rough segments. 

The scaling behaviour at the PTE critical endpoints 
in region III deserves attention. The existence of such 
points has been anticipated Q but not realized in solid- 
on-solid models before, and their scaling properties have 
not been confirmed numerically before. 

It is advantageous to align the transfer matrix along 
the preferred direction of the steps. Therefore, we switch 
our attention to the SS-rough phase boundaries where, 
inside the SS facet, the steps run parallel to the x axis, 
see Fig. [|(b). Inside this SS facet, the tilt vector is equal 
to Q = (Q x , Q y ) = (0, 1) We use the T y set-up, such that 
Q x is continuous. 

Consider the crystal shape close to the PTE point. De- 
fine q y — Q* - Q y > and q x — Q x - Q* = Q x as the 
deviation of the tilt with respect to the SS facet (Q* = 1 
and Q* = 0). Let e y = E y — E* and e x — E* — E x be 
the distances from the PTE point The crystal is rough 
(faceted) when e y > (e y < 0). The rough-to- faceted 
edge is PT-type (first-order) when e x > (e y < 0). The 
free energy f(q x ,q y ) as function of the tilt angles is an 
analytic function near PT transitions. This is a direct 
result of the dilute (free Fermion) nature of the "gas" of 
defects Jl|Jl|. Therefore we can expand it in terms of 
a polynomial in q x and q y . Our presumption is that this 
remains true also at the PTE point. 

f{q x ,q v )= ai\q x \ + he y \q y \ + 

i& 3 e.M 3 + ^ 4 M 4 + ... (5) 

Notice the absence of the quadratic term. This is a 
well known property of the entropic hardcore repulsion 
of meandering defect lines (the missing steps in the SS 
structure) at low densities [p^JT^ I. Eq. reproduces 
the usual PT scaling behaviour. The crystal shape fol- 
lows from minimizing f{q x ,q y ) with respect to q x and 
q y at constant e x > and e y > 0. This gives q x = 
and b\t y + by,e x q y + b 4 q y = 0. Along a path of con- 
stant e x > the tilt angle vanishes with a square root 
power, q y ~ |e y |2 and the free energy has the famous 
PT power law singularity, / ~ This is the power 

with which the surface orientation connects to the SS 
facet. The anisotropic scaling exponent z = 2, l y ~ I*, 
follows from the / ~ q y power with which the free en- 
ergy scales at the PT points. This implies that the 
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free energy of one single defect scales per unit length 
asm~ l x [f{2n/l x ) — /(0)] ~ l 2 x with l x a finite lattice 
length in the direction perpendicular to the defect, m is 
inversely proportional to the correlation length in the y- 
direction, m ~ Z" 1 . Our numerical results are not shown 
here, since they are in complete agreement with this. 

At the PTE point, e x = 0, the same quantities scale 
according to Eq.(||) as q y ~ | e w j ^ , / ~ q y , and z = 3. 
Along the first-order line, e x < 0, the jump in the ori- 
entation angle vanishes linearly, Aq y ~ |e x |. q x develops 
a jump as well, but with a weaker power, because the 
rotation of the surface orientation inside the rough phase 
is an higher order effect. 

It is not certain that q y is the lowest surviving power 
in the free energy at the PTE point. It could be the 
\q y \ 5 term instead. For example, in the conventional free 
fermion analysis for PT transitions the free energy is an 
odd function in | g y | . In that case the critical exponents 
change to q y ~ |e a |^, / ~ q y , and 2 = 4. 

In Fig. [l9| we show the free energy as a function of q y 
at the PTE point for K = 0.37, L = -0.26 in region III 
for system sizes 4 < N < 10. The curves are straight 
lines. This confirms that the free energy is going to zero 
as / ~ q y and therefore that z — 3. 

The scaling behaviour q y ~ across the PTE point 
and the power Aq y ~ \e x \ with which the jump in q y along 
the first-order boundary vanishes are difficult to resolve 
numerically because of the discrete nature of q y . We find 
evidence that the jump scales linearly, but the results are 
somewhat ambiguous. 

At the III to IV phase boundary the PT segment along 
the rough-SS facet boundary shrinks to zero and the two 
PTE points merge. This can be described by a polyno- 
mial of the form 

f(qx,q y ) = ai\q x \ +b 1 e y \q y \ + 

\Hc-zl)\<ly? + \H<ly\ i + --- (6) 

e x is now the distance along the rough-SS ridge from the 
point in the middle between the two PTE points, and 
c > (c < 0) in region III (IV ). This polynomial, if 
correct, would imply that the critical exponents at the 
merging point, c = 0, are the same as at the PTE points, 
z = 3 and (3 = ^, except that the jump in angle vanishes 
now as, as q y ~ | e^. | 2 . We find numerically that the free 
energy scales as / ~ q^ and, therefore, that indeed z = 3. 

VII. CONCLUSIONS 

In this paper we obtained the crystal shape for a BC- 
SOS model on a honeycomb lattice. This shape is quite 
complex. It contains two types of fully tilted facets, the 
SS and ZZ facets. The latter disappear for negative near- 
est neighbour interactions. In region III in that part of 
the phase diagram, we find that the rough-to-SS facet 
boundary becomes partially first-order with sharp edges 



where the surface orientation jumps. These first-order 
segments connect at PTE critical endpoints to conven- 
tional PT-type segments where the rough phase connects 
to the facet without a jump in orientation. The scal- 
ing properties at these PTE points follow from a simple 
higher-order polynomial free energy generalization of PT 
transitions, with anisotropic scaling with exponent z = 3 
instead of z = 2. 

In region II of the phase diagram the rough-SS phase 
boundary is PT-type everywhere, but a first-order line 
sticks out of the FRE point into the rough phase. These 
FOR lines remain very small in this model, but we re- 
solved them numerically sufficiently to be sure they exist. 

The main object of this study is to establish how gen- 
eral and universal the KPZ scaling properties of the FRE 
point in the square lattice BCSOS model are. Our con- 
clusion is that the stochastic nature of their transfer ma- 
trix makes KPZ-type FRE points special and unstable. 
In general the local crystal shape around FRE points is 
more complex than simply one facet ridge splitting into 
two PT lines. Instead there are FOR lines sticking into 
the rough phase and/or the rough- facet phase boundaries 
become first-order. Even more complex structures (not 
observed in our model) are possible as well. 

We are unable to resolve the scaling properties of the 
the critical point at the end of the FOR line (opposite 
to the FRE point). The FOR lines remain too short. 
It might still be KPZ-type scaling with z = §. In the 
exactly soluble square lattice BCSOS model, the KPZ 
nature of the scaling is closely linked to the degeneracy 
of the free energy with respect to all surface orientations 
(this is implied by the stochasticity of the transfer ma- 
trix). Our model lacks this special symmetry, but, at 
these critical points at the end of the FOR line, the free 
energies of small angle surface orientations become de- 
generate. Maybe this local version of that symmetry is 
sufficient for z = | anisotropic scaling, or it could lead to 
a new, yet unknown value of z. This aspect needs further 
study. 

This work is supported by NSF grant DMR-9700430. 



[1] G. Wulff, Z. Kristallogr. Mineral. 34, 449 (1901). 

[2] C. Herring, Phys. Rev. 82 87 (1951). 

[3] For a review see, C. Rottman and M. Wortis, Phys. Rep. 

103, 59 (1984). 
[4] C. Rottman and M. Wortis, Phys.Rev.B 29, 328 (1984). 
[5] I.M. Nolden and H. van Beijeren, Phys.Rev.B 49, 17224 

(1994) ; I.M. Nolden, J.Stat.Phys. 67, 155 (1992). 

[6] R.F.P. Grimbergen, H. Meekes, P. Bennema, 
H.J.F. Knops, and M. den Nijs, Phys.Rev. B 58, 5258 
(1998). 

[7] J. Neergaard and M. den Nijs, Phys.Rev.Lett. 74, 730 

(1995) . 



7 



[8] J.D. Shore and D.J. Bukman, Phys.Rev.Lett. 72, 604 

(1994) . 

[9] J.D. Shore and D.J. Bukman, Phys.Rev.E 51, 4196 

(1995) . 

[10] See e.g., M. Kosterlitz and D.J. Thouless, J. Phys. C 6, 
1181 (1973). 

[11] V.L. Pokrovsky and A.L. Talapov, Phys.Rev.Lett. 
42, 65 (1979); V.L. Pokrovsky and A.L. Talapov, 
Zh.Exp.Teor.Fys. 78, 269 (1980). 

[12] For a review on the Fermion approach see: M. den Nijs, 
in Phase Transitions and Critical Phenomena, Vol. 12, 
eds. C. Domb and J.L. Lebowitz (Academic Press, Lon- 
don, 1988). 

[13] M.E. Fisher, Phys.Rev.B 36, 644 (1987). 

[14] D. Dhar, Phase Transitions 9, 51 (1987). 

[15] LH. Gwa and H. Spohn, Phys. Rev. A 46, 844 (1992). 

[16] Douglas Davidson and Marcel den Nijs, in preparation. 

[17] For more detail, see e.g., M. den Nijs, J. Phys. A (1987). 



( >KPZ 




FIG. 1. Crystal structure with facet ridge endpoint as re- 
alized in the BCSOS model on a square lattice with only next 
nearest neighbour interactions. 




FIG. 2. Lattice of the triangular BCSOS model defined in 
eq.(3). 
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FIG. 3. All possible local configurations around each ver- 
tex and their Boltzmann weights. 
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FIG. 4. Phase diagram of the triangular BCSOS model at FIG. 6. Regions in the phase diagram where the various 

zero tilt, fields E x = E y = 0. crystal shapes, I-V as described in the text are realized. The 

dashed line shows the location of the roughening line in the 
zero field model. 




FIG. 7. Global crystal shape at point K — —0.24 and 
L — 0.29 in region I. 
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FIG. 10. (a) Free energy f(Qx) and (b) the vertical sur- 
face tilt Qy(Q x ) as function of horizontal tilt Q x at point 
K = -0.24, L = 0.29 in region I for E y = using the T x 
transfer matrix set-up where Q x takes only discrete values. 
Data is shown for system sizes 6 < iV < 18. 
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FIG. 11. Free energy f(E y ) as function of the surface tilt 
field E y at point K = -0.24, L = 0.29, in region I for Q x = 
using the T x transfer matrix set-up where Q y is a continuous 
variable. Data is shown for system sizes 6 < N < 18. 



1.7 



1.6 - 



N 1.5 - 



1.4 - 




0.05 0.1 0.15 0.2 

1/N 

FIG. 12. Finite size scaling behaviour of the mass gap 
scaling exponent z, at the endpoint of the FOR line, at two 
points in the phase diagram: (a) K = 1.43 and L = —0.63 
and (b) K — 0.51 and L = —0.42. The curves converge to 
values close to the KPZ dynamic exponent z = § . 
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FIG. 13. Free energy as function of the perpendicular 
tilt Q x at system size N = 18 along the FOR line at point 
K = 1.43 and L = —0.63 in region II. The value of E y is 
shown above each curve. 
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FIG. 14. The tilt angle Q y (Q x ) at system size N = 18 for 
K = 1.43 and L = -0.63 in region II when E y = 0.918. 
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FIG. 15. Finite size scaling approximations for the lo- 
cation of the endpoint of the FOR line at K — 1.43 and 
L = —0.63. The curves show the value of E y where 
Ao(Qz) = \o(Q x + V3/2N) for: (a) Q x =0 (b) Q x = V3/2N, 
and (c) Q x — 2\/3/2N. The dot shows the location of the 
FRE point. 
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g 8 1 FIG. 18. Two possible scenarios for the crossover between 

Q region II and III. 



FIG. 16. Free energy as function of the perpendicular tilt 
Q x for system size N = 18 between the FRE and PTE points 
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